function ecog = ecogModulationIndex(ecog,band, trials, chan, interval, signaltype)

% ecogModulationIndex calculates the Modulatione as proposed by Canolty_2006_Science. 
%
% INPUT: (1) ecog structure
%        (2) band: vector of numbers of specifying the band according to
%            enumeration in ecogBandpassFilter.m
%        (3) trials: vector of integers specifying the trial number desired
%        (4) interval: time of interest in seconds.
%        (5) signaltype: 'evoked' means modulation of a channel of average
%                        over all trials
%                        'trial' means modulation in each trial
% 
% OUTPUT: ecog structure with values proposed by Canolty (suppl. materials)
% 
% EXAMPLE: ecog = ecogModulationIndex(ecog,[1 6], 1:size(ecog.data,3),[8 8],[0 .5], 'trial')

clc;
switch signaltype
    case 'evoked'
        StartIdx = ecog.srate + round(ecog.srate * interval(1)) ;
        EndIdx   = ecog.srate + round(ecog.srate * interval(2)) ;
        numpoints    = length(ecog.BPdata{band(1)}(chan,StartIdx:EndIdx,1));
        skip = randsample(numpoints,(EndIdx - StartIdx - 1))+ecog.srate;

        amplitude = abs(hilbert(mean(ecog.BPdata{band(2)}(chan(2),:,:),3)));
        phase     = angle(hilbert(mean(ecog.BPdata{band(1)}(chan(1),:,:),3)));
        z         = amplitude .* exp( 1i*phase );


        ecog.m_raw     = mean(z);
        surrogate_m = zeros(length(skip),1);
        for s = 1:length(skip);
            surrogate_amp    = circshift(amplitude',skip(s))';
            surrogate_m( s ) = abs( mean( surrogate_amp .* exp( 1i*phase )));
        end

        [ecog.surrogate_mean, ecog.surrogate_std] = normfit( surrogate_m );

        ecog.m_norm_length = (abs( ecog.m_raw) - ecog.surrogate_mean)/ecog.surrogate_std;
        ecog.m_norm_phase  = angle( ecog.m_raw );
        ecog.m_norm = ecog.m_norm_length * exp( 1i*ecog.m_norm_phase );
        
        
        
    case 'trial'
        
        ecog.m_raw = zeros(1,size(ecog.data,3));
        ecog.surrogate_mean = zeros(1,size(ecog.data,3));
        ecog.surrogate_std = zeros(1,size(ecog.data,3));
        ecog.m_norm_length = zeros(1,size(ecog.data,3));
        ecog.m_norm_phase = zeros(1,size(ecog.data,3));
        ecog.m_norm = zeros(1,size(ecog.data,3));
        StartIdx = ecog.srate + round(ecog.srate * interval(1)) ;
        EndIdx   = ecog.srate + round(ecog.srate * interval(2)) ;

        for t = trials;
            numpoints    = length(ecog.BPdata{band(1)}(chan,StartIdx:EndIdx,1));
            skip = randsample(numpoints,(EndIdx - StartIdx - 1))+ecog.srate;


            amplitude = abs(hilbert(ecog.BPdata{band(2)}(chan(2),:,t)));
            phase     = angle(hilbert(ecog.BPdata{band(1)}(chan(1),:,t)));
            z         = amplitude .* exp( 1i*phase );


            ecog.m_raw     = mean(z);

            surrogate_m = zeros(length(skip),1);
            for s = 1:length(skip);
                surrogate_amp    = circshift(amplitude', skip( s ))';
                surrogate_m( s ) = abs( mean( surrogate_amp .* exp( 1i*phase )));
            end

            [ecog.surrogate_mean(t), ecog.surrogate_std(t)] = normfit( surrogate_m );

            ecog.m_norm_length(t) = (abs( ecog.m_raw) - ecog.surrogate_mean)/ecog.surrogate_std;
            ecog.m_norm_phase(t)  = angle( ecog.m_raw );
            ecog.m_norm(t) = ecog.m_norm_length(t) * exp( 1i*ecog.m_norm_phase(t) );
        end
        
end